classdef fit_simulation < handle
    %UNTITLED3 Summary of this class goes here
    %   Detailed explanation goes here
    
    properties
        i=0;
        
        %        meanThinkTime=a(1)
        %        meanCPUTime=a(2)
       mdl=[];
       % mdl_with_time=[];
       fit_=[];
       history=[];
    end
    
    methods
        function obj=fit_simulation()
            obj.mdl = @(a,x) obj.run_simulation(a(1),   a(2), x(:,1)', x(:,2)' );             
            % mdl_with_time=@(a,x) obj.run_simulation(a(1),   a(2), x(:,1)', x(:,2)' );             ;
        end
        
        function res=run_simulation(obj, meanThinkTime,   meanCPUTime, numberOfUsers, numberOfServers)
                if (meanThinkTime<0 || meanCPUTime<0)
                    res=ones(size(numberOfUsers,2),1).*10000;
                    return;
                end
                
                %[status result]=system('python ./simulator/mytest.py')
                cmdString='';
                cmdString = [cmdString, ' ',num2str(meanThinkTime)];
                cmdString = [cmdString, ' ',num2str(meanCPUTime)];
                cmdString = [cmdString, ' ',vect2str(numberOfUsers,'separator',',', 'openingDelimiter',' ', 'closingDelimiter',' ','formatstring', '%d')];
                cmdString = [cmdString, ' ',vect2str(numberOfServers,'separator',',', 'openingDelimiter',' ','closingDelimiter',' ','formatstring', '%d')]; 
                % disp(cmdString);
                cmdString = ['python ./simulator/mytest.py' cmdString];	 
                [status res_]=system(cmdString);
                % res_
                r=regexp(res_,',','split');
                res=[];                
                for i=1:length(r)
                    res=[res str2num(r{i})]; 
                end
%                 size(res)
%                 numberOfUsers
%                 numberOfServers
%                 res                
                if not(isempty(find(res>1))) 
                    res(find(res>1))=1;
                end
                    
                if (size(res,2)<size(numberOfUsers,2))
                    k=size(numberOfUsers,2)-size(res,2); 
                    res                    
                    res=[res ones(1,k).*10000];                    
                    return;
                end
                res=res';
                obj.i=obj.i+1;
                disp(sprintf('it is trying meanThinkTime=%6d and meanCPUTime=%6d at %6d', meanThinkTime, meanCPUTime, obj.i) )
                % obj.history = [obj.history res];
                % res'
        end %run_simulation       
        
        function [ahat,cov,mse]=fit(obj,a0,x,y)            
            % Fit the model to data starting from the arbitrary guess a0            
            [ahat,r,J,cov,mse] = nlinfit(x,y,obj.mdl,a0,statset('Display','iter','DerivStep', 0.1)); %eps^(1/3)))            
        end

        function x=build_input(obj,numberOfUsers,numberOfServers)
            %        numberOfUsers=x(:,1)'           
            %         numberOfServers=x(:,2)' 
            x = [numberOfUsers' numberOfServers' (1:size(numberOfUsers,2))']; 
        end
        
         function [x,y]=generate_test_data(obj)
            samples = 42;            
            numberOfUsers = round(100+25*sin(0.1*pi*1:samples)); %[52,52];
            numberOfServers= 6+round(8*rand(1,samples)); %[16,4];
            x=obj.build_input(numberOfUsers,numberOfServers);
         
            % Generate synthetic data with parameters a = [1.0; 0.10],             
            % normally distributed noise with standard deviation 0.1:
            % for reproducibility       
            rand('twister',9845);      
            MeanThinkTime = 1.0; MeanCPUTime = 0.10; 
            a = [MeanThinkTime; MeanCPUTime];
            epsn = normrnd(0,0.1,samples,1);           
            y = obj.mdl(a,x) + epsn;
            a0 = [2.0 ; 0.10];  
            [ahat,cov,mse]=obj.fit(a0, x,y)
            obj.fit_=struct('ahat',ahat,'cov',cov,'mse',mse); 
         end
        
        % clear; clear all;  se=fit_simulation; se.identifyModel
        function [x,y]=identifyModel(obj)
            load a_sample.mat
            % intrvl=(40:50);
            % intrvl=(40:89);
            intrvl=(50:218);
            sample_run_to_use=1;

            utilization=utilization(sample_run_to_use,intrvl)/100;
            servers=servers(sample_run_to_use,intrvl) ;
            num_users= num_users(intrvl); 
            % scatter3(num_users,servers,utilization); 
            x= obj.build_input(num_users,servers);
            y = utilization';
            plot(y); hold on
            % a0 = [3.0 ; 0.01];   %MeanThinkTime, MeanCPUTime = 0.10; 
            % a0 = [2.6636 0.0197];
            a0 = [6.9122    0.0158];
            [ahat,cov,mse]=obj.fit(a0, x,y)
            obj.fit_=struct('ahat',ahat,'cov',cov,'mse',mse); 
        end
    
        % clear; clear all;  se=fit_simulation;  se. show_data
        function show_data(obj)
            load a_sample.mat
            intrvl=(30:219);
            for  sample_run_to_use=1:3            
                figure;
                utilization_=utilization(sample_run_to_use,intrvl)/100;
                servers_=servers(sample_run_to_use,intrvl) ;
                num_users_= num_users(intrvl); 
                % axis([xmin xmax ymin ymax]);  
                
                set(gca,'XTick', (1:6) );
                scatter3(num_users_,servers_,utilization_); 
                set(gca,'YTick',(1:1:6))
                axis([0 900 0 6 0 1]);
            end
        end
        
        % clear; clear all;  se=fit_simulation; se.show_result(x,y) % [x y]=se.identifyModel;
        % clear; clear all;  se=fit_simulation; [x,y,yp]=se.show_result()
        % 
        function [x,y,yp]=show_result(obj)    
            load a_sample.mat
            ia=4; ib=10;
            intrvl=(ia:ib);
%           factor=1;            
%             intrvl1=(ia:(ib-ia)/(factor*(ib-ia+1)):ib); % interpolation
%             intrvl1 = intrvl1(2:end); 
            
            sample_run_to_use=1;
            utilization=utilization(sample_run_to_use,intrvl)/100;
            servers=servers(sample_run_to_use,intrvl) ;
            num_users= num_users(intrvl);             
            x= obj.build_input(num_users,servers);
            y = utilization'; 
            
           % a_hat = [2.6636 0.0197]; % i got this by solving the optimization 
            % a_hat = [6.9122    0.0158];
            % a_hat = [6.1192    0.0216]; 
            a_hat = [ 7.5068    0.0162]; 
             
            yp=obj.run_simulation(a_hat(1),   a_hat(2), x(:,1)', x(:,2)'  );                      
            % yp=obj.mdl(obj.fit_.ahat,x);
%             figure
%             scatter(y,yp); 
            figure;
            subplot(3,1,1); plot(intrvl,y,'--rs'); hold on; plot(intrvl,yp);            
            subplot(3,1,2); plot(intrvl,servers);
            subplot(3,1,3); plot(intrvl,num_users);
            
%             xrange = min(x):.01:max(x);
%             hold on
%             scatter(x,y)
%             plot(xrange,obj.mdl(obj.fit_.ahat,xrange),'r')
%             hold off
         end
    end
    
end

